# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:light
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.5
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # 任务 1.1：复势函数构造
# **依赖库声明**  
# - SymPy 1.12  
# - 默认参数：圆柱半径 `a=1.0`，来流速度 `U=1.0`

# %%
import sympy as sp
sp.init_printing(use_latex='mathjax')  # 启用 LaTeX 渲染

# 定义符号变量
x, y, a, U = sp.symbols('x y a U', real=True)
z = sp.symbols('z', complex=True)  # 独立复变量符号

# %% [markdown]
# ## 步骤 1：构造复势函数

# %%
def construct_complex_potential(U_val, a_val):
    """构造圆柱绕流复势函数
    
    参数:
        U_val (float): 来流速度
        a_val (float): 圆柱半径
    
    返回:
        sp.Expr: 复势表达式
    """
    kappa = 2 * sp.pi * U_val * a_val**2
    Phi = U_val * z + kappa / (2 * sp.pi * z)
    return Phi

Phi_symbolic = construct_complex_potential(U, a)
display(sp.Eq(sp.Symbol(r'\Phi(z)'), Phi_symbolic.simplify()))

# %% [markdown]
# ## 步骤 2：验证圆柱表面法向速度为零

# %%
# 极坐标定义
r, theta = sp.symbols('r theta', real=True)
z_polar = r * sp.exp(sp.I * theta)  # z = r e^{iθ}

# 转换为极坐标复势
Phi_polar = Phi_symbolic.subs(z, z_polar).simplify()
phi_polar = sp.re(Phi_polar)  # 速度势实部

# 计算径向速度
u_r = sp.diff(phi_polar, r).simplify()
u_r_at_a = u_r.subs(r, a).simplify()
display(sp.Eq(sp.Symbol(r'u_r|_{r=a}'), u_r_at_a))  # 输出 0

# %% [markdown]
# # 任务 1.2：解析性验证
# ## 极坐标柯西-黎曼方程验证

# %%
psi_polar = sp.im(Phi_polar).simplify()  # 流函数虚部

# 符号验证
cr1 = sp.diff(phi_polar, r) - (1/r) * sp.diff(psi_polar, theta)
cr2 = sp.diff(psi_polar, r) + (1/r) * sp.diff(phi_polar, theta)
print("极坐标 C-R 方程符号验证:")
display(sp.Eq(sp.Symbol(r'\frac{\partial \phi}{\partial r} - \frac{1}{r} \frac{\partial \psi}{\partial \theta}'), cr1.simplify()))
display(sp.Eq(sp.Symbol(r'\frac{\partial \psi}{\partial r} + \frac{1}{r} \frac{\partial \phi}{\partial \theta}'), cr2.simplify()))

# %% [markdown]
# ## 直角坐标系数值验证

# %%
phi_real = sp.re(Phi_symbolic).simplify()
psi_imag = sp.im(Phi_symbolic).simplify()

# 符号验证
cr1_cart = sp.diff(phi_real, x) - sp.diff(psi_imag, y)
cr2_cart = sp.diff(phi_real, y) + sp.diff(psi_imag, x)
assert cr1_cart.simplify() == 0, "直角坐标 C-R 方程 1 失败"
assert cr2_cart.simplify() == 0, "直角坐标 C-R 方程 2 失败"
print("直角坐标系符号验证通过！")

# %% [markdown]
# ## 数值残差验证（最大残差 <1e-6）

# %%
from utils import validate_cr_residuals

# 极坐标测试点（r > a）
test_points_polar = [
    (2.0, 0),          # θ=0, r=2a
    (1.5, sp.pi/4),    # θ=45°, r=1.5a
    (3.0, sp.pi/2)     # θ=90°, r=3a
]

# 极坐标验证（替换 a=1.0, U=1.0）
residuals_polar = validate_cr_residuals(
    phi_polar,
    psi_polar,
    variables=[r, theta],
    test_points=test_points_polar,
    extra_subs={a: 1.0, U: 1.0}  # 关键修复：替换额外符号
)
print(f"极坐标最大残差: {residuals_polar['max_error']:.2e}")

# 直角坐标测试点
test_points_cart = [
    (2.0, 0.0),    # x=2a, y=0
    (1.0, 1.0),     # x=a, y=a
    (0.5, 1.5)      # x=0.5a, y=1.5a
]

# 直角坐标验证（替换 a=1.0, U=1.0）
residuals_cart = validate_cr_residuals(
    phi_real,
    psi_imag,
    variables=[x, y],
    test_points=test_points_cart,
    extra_subs={a: 1.0, U: 1.0}  # 关键修复：替换额外符号
)
print(f"直角坐标最大残差: {residuals_cart['max_error']:.2e}")

# 断言残差 <1e-6
assert residuals_polar['max_error'] < 1e-6, f"极坐标残差超标: {residuals_polar['max_error']}"
assert residuals_cart['max_error'] < 1e-6, f"直角坐标残差超标: {residuals_cart['max_error']}"
print("\n所有数值残差验证通过！")

# %% [markdown]
# ## 渐近行为分析

# %%
V_pert = sp.diff(Phi_symbolic - U*z, z).simplify()
display(sp.Eq(sp.Symbol(r'V_{\text{pert}}(z)'), V_pert))